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A massive neutrino which decays after recombination (t > 10 13 sec) into relativistic decay 
products produces an enhanced integrated Sachs- Wolfe effect, allowing constraints to be placed on 
such neutrinos from present cosmic microwave background anisotropy data. Previous treatments 
of this problem have approximated the decay products as an additional component of the neutrino 
■ background. This approach violates energy-momentum conservation, and we show that it leads to 

Q\ ' serious errors for some neutrino masses and lifetimes. We redo this calculation more accurately, 

by correctly incorporating the spatial distribution of the decay products. For low neutrino masses 
and long lifetimes, we obtain a much smaller distortion in the CMB fluctuation spectrum than have 
previous treatments. We combine these new results with a recent set of CMB data to exclude the 
mass and lifetime range mj, > 100 eV, r > 10 12 sec. Masses as low as 30 eV are excluded for a 
narrower range in lifetime. 
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I. INTRODUCTION 



Anisotropies in the cosmic microwave background (CMB) contain an enormous amount of information about the 
universe. Data presently available has been used 0- [jll| to constrain from two up to eight cosmological parameters. 
With the promise of ever more precise measurements of these anisotropies, it has become possible to envision CMB 
fluctuations as a tool to go beyond this minimal set and constrain other areas of physics. Recent proposed constraints 
include limits on Brans-Dicke theories |Q , constraints on time- variation in the fine-structure constant JlJ|l4| , tests of 
[ finite-temperature QED , and limits on various models for both stable |l6j and unstable (T7|-[n| massive neutrinos. 
All of these additional constraints, with one exception, are based on the high-precision fluctuation spectra expected 
from the MAP and PLANCK satellites. The sole exception is reference jl^], in which Lopez et al. pointed out 
that the radiation from a neutrino decaying into relativistic decay products could produce such a large integrated 
Sachs- Wolfe (ISW) effect, that a fairly large mass-lifetime range can be ruled out from current observations. Lopez 
+3 , et al. argued that a neutrino with a mass greater than 10 eV and a lifetime between 10 13 and 10 17 sec could be ruled 
out. (Although this calculation assumes nothing about the nature of the decay products other than that they are 
relativistic, this limit is most useful when applied to decay modes into "sterile" particles such as a light neutrino and 
a Majoron, since other, more restrictive limits apply to photon- producing decays). Hannestad fl9[| showed that the 
MAP and PLANCK experiments should produce an even larger excluded region in the neutrino mass-lifetime plane. 

In this paper, we improve on a major approximation of references Jl8|,[il|. In these papers, the relativistic decay 
products were simply added to the background neutrino energy density in the program CMBFAST p0| . However, 
when the massive neutrinos decay, the spatial distribution of the decay products is determined by the distribution of 
the non-relativistic decaying particles; it is not identical to the distribution of the background massless neutrinos. In 



fact, the approach of references 1 18 19[| violates energy-momentum conservation. Although in this approach energy and 
momentum are explicity conserved at zeroth order (the mean), the first order perturbations violate energy- momentum 
conservation. This may seem like a small effect, but it actually has significant consequences for the CMB fluctuation 
spectrum. 

In the next section, we discuss the formalism for the Sachs- Wolfe effect in the presence of neutrinos decaying 
after recombination. In section III, we present our results, showing the effects of correctly incorporating the spatial 
distribution of the decay products, and provide a simple physical explanation of these effects. In section IV, we show 
how our revised calculation affects the excluded region in the neutrino mass-lifetime plane, and in section V we briefly 
summarize our conclusions. A comparison of our new results with current data leads to the excluded region rrih > 100 
eV, t > 10 12 sec, although smaller masses can also be excluded for a smaller range of r. 
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II. THE ISW EFFECT WITH AN UNSTABLE NEUTRINO: FORMALISM 



To calculate the CMB fluctuations in the presence of a decaying massive neutrino, we first review the basic precepts 
of the pertinent linear perturbation theory. The perturbed homogeneous, isotropic FRW metric can be parametrized 
as 

ds 2 = -a(r) 2 [dr 2 {l + 2ip) - dx 2 (l - 20)] , (1) 

where a is the scale factor normalized to unity today and r is the conformal time defined by dr = dt/a, t being the 
proper time of a comoving observer. This particular gauge is referred to as the conformal Newtonian gauge because 
the behavior of the potentials (</>, ip) is akin, loosely speaking, to that of the Newtonian potential. These potentials 
determine the large scale CMB behavior. In particular, the photon temperature perturbation decomposed into its 
Fourier and angular modes can be shown to be ]2lj ] 

/•TO 

A e (k) = dT(0(fc,T)+Vj(fc,T))exp(- K (T))jHfcT O -fcr), (2) 
Jo 

where the subscript '0' refers throughout to the present time and k is the optical depth from the present to some 
conformal time r in the past. For the purpose of clarity, all the sources contributing to the anisotropy from inside the 
last scattering surface have been set to zero in Eq. |[ The effect of the sources contributing to the anisotropy between 
last scattering and the present (as given in Eq. |J) is called the Integrated Sachs- Wolfe (ISW) effect. The power in the 
I th multipole is normally defined as £(£ + l)C t with @ 

C e = (4ir) 2 dkk 2 \A e {k)\ 2 . (3) 
Jo 

The sources mentioned in connection with the ISW effect can be varied. At any time, the modes that are im- 
portant to the ISW effect correspond to those scales which are smaller than the sound horizon of the whole fluid 
(matter+radiation) at that time. For these modes, the potentials can decay if there is radiation pressure or if the 
universe expands rapidly. In models with no cosmological constant, the main contribution to the ISW effect comes 
from just after recombination (since radiation redshifts faster than matter). Inclusion of a cosmological constant 
leading to a rapid expansion of the universe late in its history would boost the power on larger scales (small £). Any 
other astrophysical process which contributes to the radiation content of the universe between last scattering and the 
present will lead to an increase in the total ISW effect. One such scenario is that of a massive particle decaying around 
or after last scattering. We will consider the case of a massive neutrino decaying non-relativistically into (effectively) 
massless particles. The details of the daughter particles turn out to be irrelevant. 

To quantify the evolution of the massive neutrino density, we will consider the Boltzmann equation for its dis- 
tribution. In a homogeneous and isotropic universe, the distribution of the collisionless massive neutrino decaying 
non-relativistically into two massless particles follows |2^] 

) = _^^ ( T) (4) 

OT tdf-h 

where td is the mean lifetime of the neutrino, and eh and qh are the comoving energy and momentum: t 2 h = q 2 +m 2 a 2 , 
and a superscript '0' will be used throughout to denote unperturbed quantities. We make the following simplifications 
throughout our treatment: (1) neglect inverse decays, (2) neglect spontaneous emission, (3) neglect Pauli blocking 
factor. The solution approaches the familiar exp(— t/td) behavior as the neutrino becomes non-relativistic. The 
evolution equation for the energy density of the unstable neutrino is the integral of Eq. It reads 

pO + 3%l + P° h ) = -^<, (5) 
a td 

where overdots represent differentiation with respect to conformal time. It should be noted (as it is important if the 
decay is not completely non-relativistic) that the right hand side contains the product of rrih and n° (number density) 
and not p° h . 

We now turn on the perturbations in the metric. Although the conformal Newtonian gauge is the most useful in 
which to understand the ISW effect, for computational purposes^] we will define all our variables in the synchronous 



x The main advantage is that CMBFAST [M is written in synchronous gauge. 
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gauge. Thus, we will express the integrand in Eq. 2 in terms of perturbations in the synchronous gauge. The 
synchronous gauge has the property that the coordinate time and the proper time of a freely falling observer coincide. 
All the perturbations are in the spatial part of the metric (gtj — a 2 5ij + a 2 hij) in this gauge. The perturbation hij 
can be Fourier transformed and broken up into its trace and a traceless part as P3| 



h(k,T)kikj +6rj(k,T) ( k t kj - —8 i3 



k 2 

~ ' (6) 



Instead of working with the conjugate momentum in the perturbed space-time, we will use and eh as defined 
above pi| and in keeping with that, we will write out the perturbed massive neutrino distribution as 

f h (x,q h ,r) = fh(qh,r)[l + ^ h (x,q h ,T)] . (7) 

Due to the fact that the decay term is linear in fh, the form of the equation for the evolution of is identical to that 
of the stable massive neutrino but with /° now given by Eq. || The stable massive neutrino case has been clearly 
worked out in Ref. . 

The decay radiation rises exponentially from being negligible in the past to some maximum value at r ~ and 
then drops off as a like normal radiation. It is more informative therefore to follow the quantity r r d = P^d/ Pi> wne re 
'rd' denotes the decay radiation and pi is the cosmological density in a massless neutrino. The evolution equation for 

r rd is 

m h n° h a 

r rd = p-. 8 

The treatment of the perturbations in the decay radiation will be analogous to that of the massless neutrino as worked 
out in Ref. (23). To evolve the perturbations in the decay radiation, we will integrate out the momentum dependence 
in the distribution function by defining (in Fourier space) 

; = iMm^^^ rrd , (9) 

J dqq 3 fU d (q,T) 

where q = qn and ^ r d is defined analagously to Eq. 0. The equation governing the evolution of F r d can be worked 
out to give 

F rd + ikpF r d + 4 f | + ^±^P 2 ( M )^ r rd = r rd (n ~ MN^p,) - ^N 2 P 2 {p) + .. ^ , (10) 

!dq h qlf°( qh ,T)* h (k,c[H,T) (l - f (j^) 2 

Na(k,r) = ^ 2f0 , V '-, (11) 

Jdq h qfJ°(q h ,T) 

where p, = k ■ n and P n {p) are the Legendre polynomials of order n. The series of terms in these equations arises 
because the perturbed quantities depend on the direction of momentum and to get the contribution to a daughter 
particle with momentum q, we need to integrate over all possible q^. Thus Eq. |l0| depends on both [i and q ■ qh- 
The situation simplifies enormously for non-relativistic deca ys because each term N p: which contributes to the p th 
multipole progressively, is of 0((q^) /a p m p h ) or higher. In Eqs.jl^ and [ll], the scries has been truncated by only keeping 
terms up to 0(qf t /a 2 m 2 l ) in the i nteg rand. Similar equations for the evolution of perturbations in the decay radiation 
can be found in references [plj , |27|| Apart from iVo, the terms on the right-hand side of Eq. ( |Io| ) are completely 
negligible for non-relativistic decays. 



The use of Eq. [lOj is our only difference from the treatment in ref. |18|. In the latter paper, the relativistic decay 
products were simply added to the neutrino background in CMBFAST. This is equivalent to setting the right hand 
side of Eq. 10 to zero. Since the perturbations in the decay products are determined by the perturbations in both 
the metric and the decaying massive particles, they are correctly described by Eq. [H]. Although this may seem like a 
minor difference, it produces very large effects, as we now show. 

III. THE ISW EFFECT WITH AN UNSTABLE NEUTRINO: RESULTS 

The formalism outlined above for the evolution of an unstable neutrino and its decay products was integrated into 
the CMBFAST code p(J. We investigated a range of masses from 10 eV to 10 4 eV and lifetimes from 10 12 to 10 18 



3 



seconds. The underlying cosmology was taken to be a standard (f2 = 1) CDM model with h = 0.5 (with Hq = lOO/i 
km sec -1 Mpc -1 ); baryon density Fish 2 — .02 and scale invariant isentropic initial conditions (the same model was 
used in rcf. ]l8| ). Our results are shown in Fig. 1 for several masses and lifetimes, along with the results obtained 
by simply adding the decay products to the relativistic background. As pointed out in Ref. [^8[ there is indeed an 
enhancement in the spectra at relatively large scales due to the ISW effect produced by the decaying neutrino. We 
will see in section IV that for many values of neutrino mass and lifetime, the spectrum produced is far from that 
observed today, and therefore a large region of parameter space is ruled out due to this effect. 
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m h =10eV, t d =10 15 s 










II 

- 


_ \ 


\ 

\ 




/// 








\ 

V 




// 
// 




A 


- 










u : 














i 1 — i — i — i — h+-j 

m h =10keV, t d =10 12 s 


\ 1 


— I — I — I 


i i 1 — 


— I 1 — I — h- 


H — 










" "V 

\ 


\ 












\ 

\ 

\ 












A \ - 
• • \ 












\ 

V 














10 






100 




1000 



FIG. 1. The spectrum of CMB fluctuations for a standard CDM model (Qb = 0.08, h — 0.5) in the presence of a decaying 
neutrino with the indicated mass and lifetime (solid curve). Shown for comparison is the spectrum obtained in Ref. Jl8| where 
the decay products were added to the background neutrinos (dashed curve). The dotted curve gives the fluctuation spectrum 
in the absence of a decaying neutrino. 

The location of this ISW induced bump is determined by the lifetime of the neutrino. For lifetimes shorter than the 
age of the universe, inhomogeneities on scales k project onto angular scales I ~ kro where tq is the conformal time 
today, and we assume a flat universe. The potentials vary in time (and hence cause the ISW effect) most significantly 
at the time of decays on scales of order the sound horizon: k 2 h ~ 3/(4tJw) where w = P/p- Therefore, the bump 
in the spectrum is produced at I ~ fc s ^T ~ (ro/r,j)(4w/3) _1 / 2 . At these late times, the dominant contribution to 
w comes from the decay radiation; hence w ~ Q rd /3 where Q rc i is the fraction of critical density in decay radiation. 
Therefore, the ISW bump should be roughly at 

r d y A\l rd [T d ) 

For a matter dominated universe the conformal time and time are related as follows: r cx i 1 / 3 . For a rah — 10 eV, 
t d = 10 15 sec neutrino, fl rd — 0.15 and To/r d ~ (4 x I0 17 sec/10 15 sec) 1 / 3 ~ 7.4. Therefore, in this case we expect 
(■isw — 29. The actual peak occurs at a larger value of I, due to entropy fluctuations which decrease w, thereby 
increasing k s h and, finally, iisw- 

Notice from Figure I that we find quantitative disagreement with the results of Lopez et al. [Q (dashed curves). 
The new results show that a more accurate treatment of the spatial distribution of the decay products produces a 
surprisingly large change in the CMB fluctuation spectrum compared to the results of reference |Q . This difference 
is larger for smaller masses as can be seen in the figure. 
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At least for low masses, the most obvious difference between the the old and new spectra is the smaller size of the 
1SW effect for the new case. This difference has a physical explanation: by not properly treating the perturbations 
in the decay radiation we overestimate an important source of the potential decays that drive the 1SW effect. To 
see this we first expand the Boltzmann equation for decay radiation perturbations, Eq. ^, in multipole moments, 
Frd = Frd,i Pi, to obtain the following hierarchy, shown here for I < 2: 

S rd + l(h + 20rd) = — (S h - S rd ) , (13) 
3 V / r rd 

rd ~k 2 ( 6 -^-<j rd ) =-^9 rd) (14) 
V 4 / r rd 

cr r d - t= ( 29rd + h + 61)) = - T —a rdl (15) 
15 V / r rd 

where S rd = F rd Q /r rd , 9 rd = ikF rd i/ '4r rd and a rd — F rd ,^/2r rd . The treatment of Lopez et al. pq | is equivalent to 
neglecting the right hand sides of the equations above. This simplification breaks down near r ~ r d , where r rd /r rd is 
not negligible. 

Neglecting the r rd /r rd terms in the Boltzmann equations for the decay radiation perturbations results in errors in 
the perturbations. Let us focus on 9 rdl which turns out to be primarily responsible for the big difference. Consider 
Eq. |l4| for modes above the horizon at t ~ r d , since for these modes the approximate treatment of Ref. yl| gives 
wrong results. For these modes the k 2 terms calculated in the approximate scheme can be shown (see Appendix) to 
be roughly similar to its exact value. Then the exact solution 9 rd is related to the approximate solution 0® d by 

9rd — Qr d — 9rd ■ (16) 



where the superscript here and in what follows denotes the solution to the set of equations [13f]15| obtained in the 
approximate scheme by neglecting the feedback terms on the right hand side. The exact solution for 9 rd is therefore 
much smaller than the approximate one. Examples for several different modes are shown in Figure 0. 




FIG. 2. The difference in the variable Q r d (dipole of the decay-produced radiation) between the correct and the approximate 
treatment in Ref. Q is plotted for different modes: k~ l = 2.8 x 10 4 Mpc (solid curve), 1.4 x 10 3 Mpc (dotted curve), 73 Mpc 
(dot-dashed curve) and 23 Mpc (dashed curve). 
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These large overestimates of 9 r d lead to correspondingly large overestimates of the ISW effect and are primarily 
responsible for the differences between our spectra and those generated in Ref. fLSj . The Appendix demonstrates 
precisely how the perturbations in the decay-produced radiation affect the potentials that govern the ISW effect, and 
how treating the decay products as identical to the massless neutrinos violates energy-momentum conservation. The 
bottom line is that the ISW effect depends significantly on the behavior of 9 r d and inaccuracies in it lead directly 
to inaccuracies in the Cj's. Why does the approximation work better for higher mass neutrinos? The ISW effect is 
generated during times when the universe has appreciable radiation. For low-mass neutrinos whose decay radiation 
never dominates the energy density, the decay radiation redshifts away relative to the matter, and is only important 
near r ~ Td- Therefore, neglecting the r r d/r r d terms creates errors in the decay radiation perturbations at the crucial 
time when they are driving the ISW effect. If the neutrino is massive enough, then its decay products are important 
for a range of times with t > when the approximation is good. So the approximate treatment works better for 
higher-mass neutrinos, like the = 10 keV, td — 10 12 sec case. 

There are other visible differences between the anisotropy spectra generated in reference |Tg| ] and our more accurate 
treatment. One difference, which exacerbates the rise in power at large scales, is a drop in the small-scale ISW effect. 
For modes which enter the horizon when there is significant radiation, the Sh term in Eq. is an important source 
term. This increases S r d relative to & and since S r d is a source for the evolution of 9 r d, it implies that 9^ d < 9 r d- 
Thus there is a decrease in the ISW effect at small scales in the approximate scheme of ref . This is not visible for 
the 10 eV unstable neutrino (in Fig. |l|) because of the comparatively large signature of the first peak, but it is readily 
apparent for the 10 keV neutrino because of the large ISW effect at small scales. 



IV. COMPARISON WITH CURRENT CMB DATA 

Since the detection of anisotropies in the CMB by COBE p8[ , there have been dozens of observations of anisotropics 
on a wide variety of angular scales (refs. (2^]- [Q). We now use these observations to place more accurate limits on 
neutrino mass and lifetime. 

In ref. p8| , a very rough constraint was placed on decaying neutrino models: a model was excluded if the power 
at I = 200 was greater than at Z = 10. As we have noted in the previous section, a more accurate treatment of the 
decaying neutrinos results in a much smaller distortion in the CMB spectrum for a certain range of neutrino masses 
and lifetimes. However, as we will see, consideration of all the data leads to constraints which are almost as stringent 
as the rough contours in ref. jl8| . 

CMB experiments typically report an estimate of the band power 

where Wij is the window function which depends on beam size and chopping strategy of experiment i. Each of these 
comes with an error bar or, in the case of correlated measurements, an error matrix A/ -1 . The naive way to constrain 
parameters in a theory then is to form 

X 2 = Y1 {pi ~ C i(!°lj) M «' {pv - CV(C,)) • (18) 

2,2' 

Here we have explicitly written the dependence of Ci on the theoretical C; 's which in turn depend on the cosmological 
parameters. This naive statistic is useful only if the band power errors are Gaussian. In fact, the probability 
distribution is typically non-Gaussian, with a large tail at the high end and a sharp rise at the low end of the 
distribution. In recognition of this, and guided by some compelling theoretical arguments, Bond, Jaffe, and Knox ]4l| ] 
proposed forming an alternative statistic: 



x 2 



£ (Z, - Zi(Cif) A/#, [Z, - Z v (C,)) (19) 



where 



Zi = In {d + Xi ) (20) 
with Xi an experiment dependent quantity, determined by the noise. The covariance matrix is now 



M z = 



(pi + x^j My (pi + a^) . (21) 



G 



Bond, Jaffe, and Knox MM have tabulated and made available the relevant data from the experiments in refs. p8fl- 
pof . We use this information and formalism^ to constrain the mass and lifetime of unstable neutrinos. 

The x 2 in eq. [jl] depends on the parameters of the cosmological model. In principle, it would be nice to allow as many 
parameters as possible to vary in addition to the mass and lifetime of the neutrino. This must be balanced against the 
constraints imposed by non-negligible time needed to run the modified version of CMBFAST^. Our strategy is to vary 
the mass and lifetime of the neutrino; the overall normalization of the Cz's; the primordial spectral index (equal to one 
for Harrison-ZePdovich fluctuations); and the calibration of each experiment. For the other cosmological parameters, 
we make "conservative" choices. That is, we choose values likely to make the power on small scales (/ ~ 200) as large 
as possible compared with the power on large scales. This acts against the effect of the decaying neutrino, which 
boosts up power on large scales, and therefore leads to more conservative limits. At each point in (to, t) space, we use 
a Levenberg-Marquardt algorithm (see e.g. [^3|,Q) to find the values of normalization, spectral index, and calibration 
which minimize the x 2 defined in Eq. [l9| The contours in Fig. show these best fit x 2 in the (to, t) plane. 




A 



Log[Mass (eV)] 



FIG. 3. Contours of % 2 in neutrino mass/lifetime plane in 1 — a intervals. Inner contour is ruled out at 4<r; outermost 
contour is la. The upper right part of the figure leads to £1 > 1 (the jaggedness at logr ~ 15.5 reflects the grid size used to 
explore the parameter space). Here h = 0.5; Qb = 0.08. 

Figure || shows the constraints on the neutrino mass and lifetime for a Hubble constant h = 0.5 and Q B h 2 = .02 



2 We also account for calibration uncertainty in the manner set down in ref. |ll| . 

3 The modified version, accounting for decaying neutrinos, takes about ten times longer than the plain vanilla code. 
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in a flat (f2 = 1) matter dominated (f^A = 0) universe. The high baryon content is above the favored value of Tytler 
and Buries |0 and serves to raise the power on small scales. Masses greater than 100 eV are ruled out for almost all 
lifetimes we have explored (r > 10 12 sec). For lifetimes between 10 14 and 10 15 sec, masses as low as 30 eV are excluded 
at the two-sigma level. These results are similar to those of ref . |[8) , but more reliable because of the improvements 
in the calculated spectra and the more careful treatment of the data. 

We checked that the contours for a different set of (h, Qb) were similar to the contours in Fig|3[ Fig. |] shows the 
results for a cosmological constant-dominated universe. Again, a sizable region is ruled out, reflecting the robustness 
of the constraint. 




12 3 4 

Log[Mass (eV)] 

FIG. 4. Same as Fig. [j, but now Qa = 0.7, Ototai = 1. 

Hannestad |l9f| performed a similar calculation, using future CMB experiments to rule out decaying neutrino models, 
but he used the same approximation as in reference jl8| ; the decay products were added into the background neutrino 
density. We expect that his excluded-region contours for low masses should shrink since ISW effect is the main 
discriminator for these masses. 

It has been noted that the decay products from a very massive neutrino could keep the universe substantially 
populated with radiation or even radiation-dominated for most of its history. The presence of radiation has the effect 
of stopping the growth of density perturbations, which in a matter-dominated universe would grow as 6 ~ a. Since 
these density perturbations should (eventually) collapse into the structure we see today, it is clear that structure 
formation arguments can also provide constraints on the neutrino mass and lifetime. Very coarse constraints on the 
radiation density can placed by requiring that the scales relevant to structure formation are able to grow sufficiently 
(assuming of course, we know the initial perturbations) , as is done in Ref. p6| . In fact, for a scale-invariant initial 
spectrum, the structure formation arguments of Ref. |26| also rule out a region at the bottom-right of our excluded 



8 



region. A more detailed analysis yields more stringent constraints [£7]]. In light of this, it is important to understand 
that the constraints from CMB are most useful for low masses, i.e., for massive decaying neutrinos which do not affect 
the late-time growth of the density perturbations appreciably. Future experiments (MAP and PLANCK) have the 
potential to constrain neutrino masses as low as 1 eV and maybe even lower fl9|| . In the end, CMB and large scale 
structure constraints on massive decaying neutrinos both overlap and complement each other. 



V. CONCLUSIONS 



Our results indicate that for calculations involving the effects of decaying particles on CMB fluctuations, exact 
conservation of energy-momentum (not just conservation of the mean energy-momentum) is crucial. When pertur- 
bations in the decay products are correctly treated as being determined by the perturbations in both the metric and 
the decaying massive particle, energy and momentum of the massive particle plus its decay products are conserved. 
The result is a much smaller change (when an unstable neutrino is added) in the CMB fluctuation spectrum than was 
noted in ref. fig] . However, by using a comparison with current data, rather than a simple constraint on C200/C10, we 
have been able to obtain an excluded region only slightly less restrictive than that obtained in ref. |l8| . This excluded 
region will grow as more data becomes available, culminating potentially in very restrictive limits from MAP and 
PLANCK pj| . Our results, of course, can be generalized to arbitrary decaying particles. 

The CMB spectra used in this work were generated with a modified version of CMBFAST pc|| . We thank Lloyd 
Knox for providing the data used to generate the constraints in section IV. This work was supported by the DOE 
and the NASA grant NAG 5-7092 at Fermilab and by the DOE grant DE-FG02-91ER40690 at Ohio State. 



APPENDIX A: SOURCE OF DECAYING POTENTIALS 



Here we show that the ISW effect in the decaying neutrino model is primarily driven by the dipole of the decay- 
produced radiation, 9 r d- The ISW effect is driven by time changes to the potentials (eq. 2), which are determined 
from Einstein's equations. In synchronous gauge, the source of these time changes is 



■ + tjj = ti + t 2 + t 3 + t 4 



(Al) 



where 



ta 



-2-77, 
a 



U = 



2fc 2 



(5 + 3w)h, 



fir V a 



(A2) 



The quantity D a is related to the anisotropic stress of the fluid: D a = —(3/2) (a/a) 2 (1 + w) a, where w = P/p 
is the equation of state of the universe. We will consider the behavior of superhorizon-scale perturbations, where 
kr -C L We assume that the neutrinos decay well into the matter-dominated phase of the universe, and that the 
decay radiation never dominates the energy density of the universe, but does come to dominate the standard radiation, 
i.e., photons and massless neutrinos. Then the equation of state takes a simple form near neutrino decay: w ~ 1/3 fl r d- 
In addition, the total fluid perturbation sources 9 and a are dominated by the decay radiation, so that we can write 
9 ~ 4w9 r d, and a ~ 4w<r r d- These assumptions are well motivated for m/j = 10 eV, td — 10 15 sec neutrinos which 
decay well into the matter dominated era, with fl r d — 0.15 at decay. 

We first examine the behavior of (fi + ip in the approximation where we neglect the r r d/i'rd terms in the Boltzmann 
equations for the decay radiation perturbations, i.e., we treat the decay radiation as massless neutrinos as in Lopez, 
et al ]T^ ] . We then consider the effect of relaxing the approximation and calculating the decay radiation perturbations 
correctly. We denote the use of the Lopez et al. []l8|| approximation in all quantities by the superscript-a. 

The potentials do not decay in a completely matter-dominated universe; <fi + "0 is sourced by the decay radiation 
and is therefore first order in w. The term t\ is directly related to a and so is of order w. The linearized Einstein 
equations imply that f) oc 8 oc w 9 r d, so that t\ is also of order w. However, t% and if are each zeroth order in w, so 
their sum must cancel to lowest order. Using the linearized Einstein equations and the continuity equation we find 
that 



(A3) 







demonstrating the required cancellation. In our approximation, the decay radiation perturbations can be calculated 
from the Boltzmann equation for massless neutrinos, which admit analytic solutions for the superhorizon modes of 
9^ d and a^ d . Using these solutions, it can be shown that 

^-29^,^12^, (A4) 

T T 

so that i", + £3 and t\ each contribute roughly comparable amounts to 4> + ip . In calculating the effect of the 
approximation on 4> + ip we will therefore have to separately consider each term. The quantities D a and 77 are very 
much affected by the approximation, since they directly depend on the decay radiation perturbations, and the error 
in the decay radiation perturbations is of order the quantities themselves. This implies that St) ~ \fj\ and 5D a ~ \D a \, 
where 8x = \x — x a \ is the absolute error in the variable x. The zeroth order quantities 77 and h are much less affected 
by the approximation. For super-horizon modes, we do not expect 77 to evolve much from its initial value, and so the 
error in it (determined by the error in 77) is naturally small. Following this line of reasoning, one can write 

Si] ~ <5?7 t ~ k 2 t 2 w 77, (A5) 

Sh ~ fc 2 t St) ~ fc 4 r 3 w 77. (A6) 

Using these relations we find that 

St!-^, St* ~ (fcr) 2 ^ , 8t 3 ~ (fcr) 2 ^ , ^ ~ O , (A7) 

T T T T 

which makes it clear that for superhorizon modes, the errors in t\ and t\ dominate the error in § + if> . Numerically 
it is seen that the error in t\ is the most important. 

We can see why the error in t\ is the most important in a simple way. The dominant source term (see eq. [IjJ ) for 
®rd i s Kd- N° w S r d and 5^ d differ only by the terms on the right-hand side of eq. [l3] (since h is not much affected by 
the approximation and 9 rd -C h for super- horizon modes). But the right-hand side of eq. [l^ contains the difference 
of Sh and S rd , and hence the fractional error in S rd is expected to be much smaller relative to that in 9 rd or a rd . So 
we can assume that S r d and 5® d are roughly the same for the purpose of estimating the errors in 6 rd and a r d in the 
approximate scheme. Therefore, for super- horizon modes (at r ~ r d ), we can write to a good approximation 

9 rd ~ ¥ rd - — 9 rd and a rd ~ &? d - — a rd . (A8) 

T*rd Trd 

From this equation, we can gauge that the fractional errors in 8 rd and a rc i are ro ughl y the same, and close to 



-Tr r d/r r d- But since the coefficient for ti in (f> + tp is much less than that for t± (see eq. A4), we expect that the error 



in t\ dominates which implies that 8{<f> + t/>) cx 77 oc 9 rd . From eq. A8 we have that > 9 rd for super- horizon modes 
at t ~ r d . Therefore for these modes, 



(j) a + yj a > <j> + ^ 



(A9) 



This is the reason for the dramatic rise in power at large scales when we neglect the decay terms in the decay radiation 
Boltzmann equations. 

Merely adding the decay radiation to the massless neutrino background causes large errors in the ISW effect. 
However, there exists a method of calculation that yields good results without introducing a separate Boltzmann 
hierarchy for the decay radiation. Since this method might be helpful for other late time processes which affect the 
CMB spectrum, and since it demonstrates that we really have isolated the source of our disagreement with Lopez et 
al. [Q, we present it here. 

The fix can be accomplished by explicitly evolving a = l/(2k 2 )(h + 677), a quantity that contains the problematic 
77, within CMBFAST in the following differential equation: 

a + 2-a = T7- (-') (l+w)a (A10) 

The dominant quantity on the right-hand side is 77, which is quite unaffected by the approximation (recall that 
5i] ~ fc 2 t 2 wq. In contrast, when a is set by the equation (the default in CMBFAST), 

a 3 a „ 9 / a x 
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the 6 term is important, and inaccuracies in the l arge scale behavior are generated. The condition that a given by 



Eq. ( A10 ) should match that obtained from Eq. ( All ) is conservation of momentum for the massive neutrino and 
its decay products. In the approximate scheme, the unperturbed quantities for the decay radiation are calculated 
correctly while the perturbations in it are set equal to that of the massless neutrino. This violates the energy- 
momentum conservation conditions for the system of the massive neutrino plus its decay products, and this is the 
reason behind the fact that different combinations of the Einstein equations lead to different potential decay rates. 
It may be noted that the CMBFAST code |M| used to calculate the fluctuation spectrum implicitly assumes energy- 
momentum conservation. When this condition is violated, the code cannot produce internally consistent results. 
In Fig. ||, we have plotted the result of using Eq. A10 (in place of Eq. All) with the approximate scheme and as 
expected, good agreement with the actual curves is obtained. This exercise clearly shows that it is important to check 
for energy-momentum conservation when using approximate methods to model any part of the energy-momentum 
tensor. 




FIG. 5. The spectrum of CMB fluctuations for a standard CDM model (fif, = 0.08, h — 0.5) in the presence of a decaying 
neutrino with the indicated mass and lifetime (solid curve). The dashed curve is the spectrum obtained when the decay 
products are added to the background neutrinos and a is evolved as an explicit differential equation instead of being fixed as 
in the standard CMBFAST code. The dotted curve gives the fluctuation spectrum in the absence of a decaying neutrino. 
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